Abelian Manna model in three dimensions and below 
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The Abelian Manna model of self-organized criticality is studied on various three-dimensional 
and fractal lattices. The exponents for avalanche size, duration and area distribution of the model 
are obtained by using a high-accuracy moment analysis. Together with earlier results on lower- 
dimensional lattices, the present results reinforce the notion of universality below the upper critical 
dimension and allow us to determine the coefficients of an e-expansion. By rescaling the critical 
exponents by the lattice dimension and incorporating the random walker dimension , a remarkable 
relation is observed, satisfied by both regular and fractal lattices. 
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INTRODUCTION 



Critical phenomena play an important role in our un- 
derstanding of complex systems in nature. One of the 
key features of critical systems is the notion of univer- 
sality It suggests common underlying mechanisms 
of seemingly different phenomena and also allows us to 
study complicated natural systems through analysis of 
simple (numerical) models. A huge number of numerical 
models have been proposed to study different features of 
critical systems. Traditionally, those systems require ex- 
ternal fine tuning of control parameter to critical point . 
However, many models exhibit the features of the critit- 
cal state without the need of tuning a control parameter, 
which is known as self-organized criticality [2j. While 
models of traditional critical phenomena have been very 
well studied both analytically and numerically, and im- 
portant results including exact ones have been obtained 
[e-<?.0], little success has been achieved for self-organized 
critical phenomena. One example of such models is the 
Manna model J|| which has so far defied any attempt for 
an analytical approach, but has been studied extensively 
numerically. Yet the question of universality seems to 
have been overlooked in the literature in a very funda- 
mental point: Does the Manna model display the same 
critical behaviour on different lattices of the same di- 
mension? Recent extensive numerical studies (BJ of the 
Abelian version of the Manna model JH [H, 0] on vari- 
ous lattices in integer dimensions 1 and 2 provide very 
strong support for the model's universal behaviour. The 
situation, however, is more complicated in non-integer 
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dimensions [8|. Thus, it would be of great interest to see 
if the results in fractal dimensions can be reconciled with 
those in integer dimensions in a systematic manner. 

Our motivation for this study is to provide a complete 
numerical picture of the Manna model. Three main re- 
sults are reported: Firstly, we confirm universality in 
three dimensions, i.e. the independence of critical ex- 
ponents and moment ratios from the detailed structure 
of the underlying lattice. This allows us, secondly, to 
firmly estimate the coefficients of an e-expansion for the 
exponents. Thirdly, we identify a general scaling relation 
unifying critical behaviour on regular and fractal lattices. 



2. MODEL AND OBSERVABLES 

The Abelian Manna model [H-0] is defined on a lattice 
of dimension d with iV sites and linear extent L, where 
N oc L d asymptotically. Each site i has a local degree of 
freedom Zi > which can be thought of as the number 
of particles residing at that site. If Zi > 1, the site is 
said to be active or unstable, otherwise it is stable. The 
system evolves by driving and subsequently fully relax- 
ing it. Driving: The system is "charged" by picking a 
site j randomly and uniformly and increasing Zj by 1. 
Relaxation: An unstable site i is picked randomly and 
uniformly; its particle number is reduced by 2 and the 
particle number of two of its qi neighbours, which are 
chosen independently and at random (possibly the same 
one twice), is increased by 1, thereby possibly render- 
ing them unstable. Dissipation takes place only when 
relaxing boundary sites transfer particles to q( v ' virtual 
neighbours "outside" the lattice. These virtual neigh- 
bours cannot topple themselves and are chosen so that 
the topology of the finite lattice corresponds to that of an 
"offcut" from an infinite lattice. Relaxation of unstable 
sites continues until there are no unstable sites left. Only 
then the system is driven again, known as a separation of 
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time scales. The number of topplings between two driv- 
ing steps is the size s of the "avalanche" and the number 
of distinct sites toppling is the area a. The duration t of 
the avalanche is measured on the microscopic time scale, 
which advances in steps of 1/N a , where N a is the instan- 
taneous number of active sites, mimicking a Poissonian 
decay of active sites. The moments of the observables 
mentioned above are measured in the stationary state, 
which is reached after a generously estimated transient. 
A number of (asymptotic) key characteristics of the lat- 
tices are listed in Table [JJ such as the average number of 
neighbours q, the average number of virtual neighbours 
g( u ) among sites with at least one virtual neighbour, and 
the particle density £ in the stationary state. 

The probability densities V^ x ' (x, L) of the observables, 
x G {.s, a, t}, are expected to display finite size scaling 

V M{x,L)=a x x-T*Q x (^L^, (1) 

provided that L>1 and x ^> x for some lower cutoff 
xq. The metric factors a x and b x are not expected to 
be universal H, whereas the exponent D x should only 
depend on the dimension and t x on the dimension and the 
boundary conditions fioj . The universal scaling function 
Q x is characterised below by moment ratios. For n > 
t x — 1, the moments scale asymptotically like (x n ) cx ZA> 
with nffl = D x (l + n — t x ). For historic reason, the 
exponents are denoted as D (for D s ), t (for r s ), z (for 
D t ), and a (for r t ). 

Five three-dimensional and five fractal lattices are em- 
ployed in this study. The three-dimensional lattices [ll| 
arc built upon the standard simple cubic (SC) lattice. 
The body centered cubic (BCC) and face centered cubic 
(FCC) lattices are also studied with next nearest neigh- 
bour interactions (BCCN and FCCN, respectively). The 
total number of sites N of all five lattices are chosen to 
be as close as possible to one another. Typically, six sys- 
tem sizes ranging from N = 181 3 to N = 1024 3 are used. 
A total of approximately 50,000 CPU hours have been 
spent on three-dimensional systems. 

The key features of the fractal lattices are listed in 
Table [Til Of those the lesser known semi-inverse square 
triadic Koch (SSTK) lattice [HEl has Hausdorff dimen- 
sion d = In 5/ In 3 and is exemplified in Fig. [TJ Of partic- 
ular interest is the Sierpinski tetrahedron (SITE) lattice, 
which is the three-dimensional version of the well-known 
Sierpinski gasket, based on tctrahedra instead of trian- 
gles. Its fractal dimension is d = 2 and thus allows a di- 
rect comparison to regular two-dimensional lattices. The 
strongly anisotropic extended Sierpinski gasket (EXGA, 
d= In 6/ In 2) is obtained by stacking L copies of a Sier- 
pinski gasket on top of one another and applying periodic 
boundary conditions. Finally, the arrowhead (ARRO, 
d = In 3/ In 2) and the crab (CRAB, d = In 3/ In 2) lat- 
tices are the same as the ones used in Q . Typically, four 
system sizes corresponding to iterations from 6 to 9 are 
used for all fractal lattices (5 to 8 for SSTK). A total 
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FIG. 1: Semi-inverse square triadic Koch lattice at iteration 
3. 

of approximately 75,000 CPU hours have been spent on 
fractal lattices. 



3. RESULTS 

Details of the Monte Carlo simulation, the fitting pro- 
cedures and the derivation of the error bars can be found 
in [BJ. In short, individual moments (x n ) (L) are fit- 
ted against the system size N oc L d to obtain the scal- 
ing exponent /Ltn /d using a power law with corrections. 
For regular three-dimensional lattices, the form of fitting 
function used for avalanche size and duration is 

{X n ) (AO = A^N^/d + B (x) N ^/d-0.25 ( 2 ) 

and for avalanche area is 

(a n )(N) = A^N^/'t + B^N^/' 1 - - 25 

+ C { n a) N^ a)/d -°- 5 . (3) 

For fractal lattices we use Eq. (2) in Q 

(x n ){N) = A ( £ ) N^ ) ' d + B<£ ) N^ )/d - 1 

+ C^N^' d - 2 . (4) 

( x ) 

The estimated scaling exponents ^„ are then linearly 
fitted against the moment orders n (n = 2,3,4,5 for all 
except for /i^ in three dimensions with n = 2, 3, 4). The 
slope gives the exponent D x and the interception with 
the abscissa gives the exponent t x , except for t s of three- 
dimensional lattices, whose estimate is obtained by em- 
ploying the exact relation D s (2 — t s ) = 2. Although the 
relative errors are as small as 3 x 10 -3 , we were unable 

(s) 

to adjust the fitting scheme as to recover fi\ = 2 within 
less than 3 standard deviations. 

The quality of all data fitting reported in this work 
is assessed by the goodness-of-fit q [14| , which is consid- 
ered good if q > 0.1, otherwise they are marked by [•]. 
Tables UlTI and FVl summarise the estimated critical expo- 
nents, all obtained with q > 0.5. For regular lattices, 

our results compare well with the literature muni (also 

[l8l Qjl for absorbing state phase transitions) , although 
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TABLE I: Key features of the three-dimensional (d = 3) lat- 
tices considered in this work. The random walker dimension 
d w is 2 for every regular lattice [2(j. The asymptotic site av- 
eraged number of nearest neighbours is q, with boundary sites 
having on average q( v ) neighbours missing (virtual neighbours 
The stationary particle density, equal to the average den- 
sity of (singly) occupied sites, is given by £. 



Lattice d d w q q( v ) ( 
SC 3 2 6 1 [0.622325(1)] 
BCC 3 2 8 4 [0.600620(2)] 
BCCN 3 2 14 5 [0.581502(1)] 
FCC 3 2 12 4 [0.589187(3)] 
FCCN 3 2 18 5 [0.566307(3)] 



some variability and discrepancy is observed in particular 
for z which may be explained by the use of slightly differ- 
ent model definitions (and dynamics) by other authors. 
A number of scaling relations (see below) are confirmed, 
such as -E x = D x (t x - 1) = B y {r v - 1) = -E y [11). 
Overall estimates are included in Table IIII1 the correla- 
tion of E s , S t and S a is taken into account by multiplying 
their respective errors by \/3- 

The moment ratios 0, [2l[ g n x ^ (to leading order) are 
independent of the system size and characterise the scal- 
ing function Q x . The fitting of moment ratios follows a 
similar procedure as the avalanche exponents, using 
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Ai x) +Bi x) N-°- 25 +Ci x} N-°- 5 . 



(5) 



Together with the avalanche exponents, they provide 
very strong support for universality in regular lattices. 
Table IIVI lists the overall moment ratios based on five 
three-dimensional lattices. 

Surprisingly (see Q for the same phenomenon in one 
and two dimensions), the amplitudes of the leading order 
of the moments of the avalanche area seem to be univer- 
sal. We found 



(a 1 ) = [0.202(4)]7V« 

(a 2 ) = 0.0151(15)^2°' 

(a 3 ) = 0.0027(6)iV' i 3 a) 

(a 4 ) = 0.00055(19)^° 

(a 5 ) = 0.00012(6)^5°' 



(6) 



with universal fi n = n + 1 — 1.4396(8) across the three- 
dimensional lattices introduced above. It is obviously 
crucial to consider (a") as a function of N, as fitting 
against L = XN 1 ^ leads to different amplitudes, because 
A varies from lattice to lattice. 



3.1. Regular lattices 

All critical exponents including previous results Q are 
summarised in Table IVIl Firstly, on regular lattices, a 



TABLE II: Key features of the fractal lattices studied, as 
listed in Table U for the three-dimensional ones. The random 
walker dimensions are exactly known or derived (in case of 
SSTK estimated) using the methods described in 



Lattice d 



d v 



q qW C 



SSTK 1.464. 
ARRO 1.584. 
CRAB 1.584. 
SITE 2 
EXGA 2.584. 



2.552. 
2.322. 
2.578. 
2.584. 
2.321. 



3 1 
7/3 1 
7/3 1 
6 3 
6 2 



[0.8435(2)] 

[0.862(2)] 

[0.8794(6)] 

[0.7427(3)] 

[0.65640(8)] 



relation between D x , t x and the dimension d can be ob- 
tained by fitting exponents against a proposed function 
D x = f x (d) and r x = h x (d). With six exponents six func- 
tions are to be determined, which, however, are related 
by scaling laws. They are 



D(2-r) = 2 



(7) 



on regular lattices (exact 



, D a = d (assumed to hold 
on regular lattices by fl5, |24|. and in the present case 
confirmed for fractal lattices) and D x (t x — 1) = — Yi x with 
Tj a = Yj s = Tj t (narrow distribution assumption f25j). 
Using r = 2 - 2/D, D a = d, r a = 1 + (D - 2)/d, and 
a = 1 + (D — 2)/z, there are thus only two functions to 
determine, which are best expressed in terms of e = 4 — d 
since d c = 4 is the upper critical dimension where 



the exponents are known exactly. Writing D = 4 - 



(s) 

., at most two amplitudes c\ can reasonably be 
determined on the basis of the three data points available. 
A fit of D with only a linear term produces a very poor 
goodness-of-fit, which does not improve satisfactorily by 
including a term quadratic in e. Omitting the quadratic 
gives 



D = 4 - 0.654(6)e + 0.0079(10)e 3 



(8) 



q « 0.095 (c[ s} = 0.60(4), 
= —0.019(7) with three terms) 



with 

— — u.uxi7 V i; vviuii 

2-0.239(4)e+0.0056(6)e 3 
goodness-of-fit. 



2 - -0.05(3), 

Similarly, z = 
however with nearly vanishing 



3.2. Fractal lattices 

Attempting to unify the above e-expansion obtained 
for regular lattices with the results for fractals with Haus- 
dorff dimension d is bound to fail, which is immediately 
clear when comparing the exponents found for the fractal 
SITE lattice (d = 2, Table E} with those for the regular 
two-dimensional lattices, or the ARRO and the CRAB 
lattice, Table [V] which have the same Hausdorff dimen- 
sion. As is well understood, the basic scaling relation 
Eq. (I?]) is valid only for regular lattices and has to be 
generalised to 



D(2 - r) =d n 



(9) 
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TABLE III: Avalanche exponents of five three-dimensional lattices. The estimates for r and D(t — 1) are obtained from D via 
the exact scaling relation D{2 — r) = 2. Identities D a — d and jj,^ = 2 are used to validate the fitting scheme. 



Lattice 


d d w 


£> 


t z a D a 


(s) 

Ta Ml 


— S s —St — S a 


SC 


3 2 


3.38(2) 


1.408(3) 1.779(7) 1.784(9) 3.04(5) 


1.45(4) 2.0057(5) 


1.38(2) 1.395(16) 1.36(13) 


BCC 


3 2 


3.36(2) 


1.404(4) 1.777(8) 1.78(1) 2.99(2) 


1.444(18) 2.0030(5) 


1.36(2) 1.390(19) 1.33(6) 


BCCN 


3 2 


3.38(3) 


1.408(4) 1.776(9) 1.783(11) 3.01(3) 


1.44(3) 2.0041(6) 


1.38(3) 1.39(2) 1.32(7) 


FCC 


3 2 




1.402(8) 1.765(16) 1.78(2) 3.1(2) 


1.48(14) 2.0035(11) 




FCCN 


3 2 


3.38(4) 


1.408(7) 1.781(14) 1.787(18) 3.00(4) 


1.44(3) 2.0051(8) 


1.38(4) 1.40(3) 1.32(9) 


Overall 3 2 


3.370(11) 


1.407(2) 1.777(4) 1.783(5) 3.003(14) 1.442(12) 2.0042(3) 


1.380(13) 


[15] 


3 2 


3.33 


1.43 1.8 






[16] 


3 2 


3.302(10) 


1.713(10) 






[17] 


3 2 


3.36(1) 


1.41(1) 1.76(1) 1.78(2) 






[22] 


3 2 




1.41(2) 1.823(23) 1.77(4) 







TABLE IV: Overall estimates of moment ratios for three- 
dimensional lattices. 



Observable x g% 



(x) 



(x) 

91 



» 



Ix) 
96 



Size s 2.373(16) 7.76(17) 30.0(14) 121(8) 

Duration t [4.164(6)] [25.99(9)] [201.4(12)] 1811(18) 
Area a 2.331(4) 7.30(5) 27.1(3) 113(2) 



with random walker dimension d w > 2 [8j, |2C§ . 

It turns out, however, that D is essentially a linear 
function of d and d„, 



D = ad + bd n 



(10) 



with the same coefficients a and b for both, regular 
and fractal lattices, which can be extracted from the e- 
cxpansion obtained above with d c — 4 = 2d w , because 
d w = 2 on regular lattices, a = c± and b — 2(1 
so that on the basis of Eq. 
0.692(12) or a = 0.60(4) and b 



l 

a = 0.654(6) and b = 
= 0.80(8) depending 



Fitting the data in Table [VI] against Eq. ([T 



on c\' . 

gives a = [0.550(4)], b = [0.822(3)], and a fitting with 
the constraint b = 2(1 — a) from exactly known val- 
ues of exponents at d = d c = 4 gives a = [0.6061(5)], 
b = [0.7878(10)]. Figure [U shows D/d as a function of 
d w /d for all lattices listed in Table [VTl As expected from 
Eq. ([T0|) . fractal and regular lattices display essentially 
the same linear dependence D/d = a + bd w /d. Above 
the upper critical dimension d c = 4 regular lattices at- 
tain their mean- field values, D = 4, r = 3/2 [23| and 
d w = 2, therefore following Eq. (fl"0|) with a = and b = 2. 
One may wonder whether fractal lattices with Hausdorff 
dimension greater than d c = 4 have correspondingly ex- 
ponents D = 2d w . 

The choice of rescaling exponent D by dimension d of 
the lattice is not random, but rather a natural choice, 

(x) 

given that we perfomed all fitting of ^„ against the 
number of sites rather than the lattices' linear length 
L and multiplied the results by the Hausdorff dimension 
d. The gap exponents for the scaling in N is D jd, which, 
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FIG. 2: The data of Table IVTl plotted in the form D/d versus 
dw/d as suggested by Eq. (|10[l . The dashed straight line is 
based on the estimates a = [0.6061(5)] and b = [0.7878(10)], 
the dotted line is the mean field theory, a = 0, b — 2. 



as it turns out, displays a very systematic dependence on 
dw/d. 

Further investigation shows that D / d fits very well to 



with a = 1.020(2) and b 
results in 



(r-5) 



(11) 



0.481(3) for all lattices which 



D = 4— 0.658(5)e + 0.00962(13)e 2 + 0.00161(3)e 3 + • • • 

(12) 

using D(2 — r) = d w = 2 for the regular ones. 

The form of Eq. (fTTj) was obtained by first fitting r 
against (D/d) K , which gives a k deviating from —2 by less 
than 2%. The coefficient a and b are then fitted according 
to Eq. (fTTj) . Figure [3] compares that relation to results for 
lattices in all dimensions. In the same manner, a similar 
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TABLE V: Avalanche exponents of five fractal lattices. 



Lattice d 



flu 



D 



w 
Mi 



SSTK 1.464.. 2.552.. 2.94(3) 1.13(2) 1.817(17) 1.21(2) 1.466(5) 1.273(11) 2.551(6) 0.37(6) 0.38(4) 0.399(17) 

ARRO 1.584.. 2.322.. 2.793(2) 1.173(2) 1.673(1) 1.280(2) 1.5847(3) 1.2985(6) 2.3103(4) 0.484(5) 0.468(3) 0.473(1) 

CRAB 1.584.. 2.578.. 3.020(5) 1.151(4) 1.837(3) 1.237(4) 1.5847(8) 1.279(2) 2.5655(12) 0.456(11) 0.435(7) 0.443(3) 

SITE 2 2.584.. 3.232(6) 1.211(4) 1.870(4) 1.357(4) 1.9975(9) 1.339(2) 2.5533(6) 0.682(14) 0.667(8) 0.677(3) 

EXGA 2.584.. 2.321.. 3.352(4) 1.312(3) 1.835(3) 1.581(3) 2.5895(6) 1.3915(8) 2.3000(2) 1.046(10) 1.066(6) 1.014(2) 



TABLE VI: Summary of exponents in all dimensions. 



Lattice d 


dw 


D 


T 


2 


a 


D a 




0) 


— E 


Ref. 


regular 


1 


2 


2.253(14) 


1.112(6) 


1.445(10) 


1.18(2) 


0.998(3) 


1.259(11) 


1.996(3) 


0.26(2) 


® 


SSTK 


1.464. 


. a 2.552.. 


2.94(3) 


1.13(2) 


1.817(17) 


1.21(2) 


1.466(5) 


1.273(11) 


2.551(6) 


0.40(3) 


here 


ARRO 


1.584. 


. b 2.322.. 


2.7938(19) 


1.1731(16) 1.6732(12) 


1.2797(17) 


1.5847(3) 


1.2985(6) 


2.3103(4) 


0.4730(16) 


hcrc+[8] 


CRAB 


1.584. 


. b 2.578.. 


3.020(5) 


1.151(4) 


1.837(3) 


1.237(4) 


1.5847(8) 


1.2793(17) 


2.5655(12) 


0.442(4) 


here+[8j 


regular 


2 


2 


2.750(6) 


1.273(2) 


1.532(8) 


1.4896(96) 


1.995(3) 


1.382(3) 


1.9993(5) 


0.761(13) 


[5] 


SITE 


2 C 


2.584.. 


3.232(6) 


1.211(4) 


1.870(4) 


1.357(4) 


1.9975(9) 


1.3388(14) 


2.5533(6) 


0.676(5) 




EXGA 


2.584. 


. d 2.321.. 


3.352(4) 


1.312(3) 


1.835(3) 


1.581(3) 


2.5895(6) 


1.3915(8) 


2.3000(2) 


1.020(3) 




regular 


3 


2 


3.370(11) 


1.407(2) 


1.777(4) 


1.783(5) 


3.003(14) 


1.442(12) 


2.0043(3) 


1.380(13) 


here 


regular 


4 e 


2 


4 


1.5 


2 


2 


4 


1.5 


2 


2 


[23] 



a hi5/ln3. b ln3/ln2. c Fractal lattice. d ln6/ln2, strongly anisotropic. 

e At upper critical dimension d c — 4 [191 . exponents take mean-field value [23]. 



relation can be obtained for z and a, 

-a) =6 (13) 

with [a = 0.936(2)] and [b = 0.3768(12)]. 

The above results suggest that the scaling in N is more 
suitable for fractals than the scaling in L. We suspect 
this is related to L not capturing the chemical distance, 
which is the distance particles need to travel on the lat- 
tice, whereas L is measured as a Euclidean distance. By 
using d w , which is sensitive to the chemical distance, and 
considering the scaling against N, which is a well-defined 
measure of the size for any lattice, we are able to deter- 
mine the relations above. 



aim to provide a complete picture about the model be- 
low the upper critical dimension. The results confirm the 
consistent and robust universal behaviour of the Manna 
model across different, regular lattices, which allowed us 
to produce an e-cxpansion of avalanche exponents below 
the upper critical dimension. A relation between crit- 
ical exponents and lattice dimension is observed which 
systematically reconciles integer dimensional with frac- 
tal lattices. 



4. CONCLUSION 

In conclusion, we studied Abelian Manna model on 
various three-dimensional and fractal lattices with the 
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